## Load libraries


suppressPackageStartupMessages({
library(ArchR)
library(rhdf5)
library(tidyverse)
library(reticulate)
library(zellkonverter)
library(Matrix)
library(dichromat)
#library(caret)
h5disableFileLocking()})
proj <- loadArchRProject("15_gene_scores_from_p2g_as_gene_expr_matrix/")

marker_genes <- c("Lamb1",  "Sparc", "Elf5", "Ascl2", "Tfap2c", "Ttr",
                  "Apoa2", "Apoe", "Cystm1", "Emb", "Spink1",  "Krt19",
                  "Dkk1", "Grhl3", "Trp63", "Grhl2",  "Pax6", "Pax2",
                  "En1", "Foxd3", "Tfap2a", "Pax3", "Sox9",
                  "Six3", "Hesx1", "Irx3", "Sox2", "Hoxb9", "Cdx4", 
                  "Hes3", "Hba-a2", "Hba-a1",  "Hbb-bh1", "Gata1", "Cited4",
                   "Cdh5", "Pecam1", "Anxa5", "Etv2", "Igf2",
                  "Krt8", "Krt18", "Pmp22", "Ahnak", "Bmp4", "Tbx4", "Hoxa11", 
                  "Hoxa10", "Tnnt2", "Myl4",  "Myl7", "Acta2", 
                  "Smarcd3", "Tcf21", "Tbx6", "Dll1", "Aldh1a2", "Tcf15", 
                  "Meox1", "Tbx1", "Gbx2", "Cdx1", "Hoxb1", "Hes7", "Osr1", 
                  "Mesp2", "Lefty2", "Mesp1", "Cer1",  "Chrd", "T", 
                  "Foxa2", "Pax7", "Fgf8", "Lhx1", "Gsc", "Mixl1", "Otx2", "Hhex",
                   "Ifitm3", "Nkx1-2", "Eomes", "Nanog", "Utf1", 
                  "Epcam", "Pou5f1")

Gene Scores UMAPs


#marker_genes <- c("Hba-a2", "Hba-a1",  "Hbb-bh1")
p <- plotEmbedding(
    ArchRProj = proj, 
    colorBy = "GeneExpressionMatrix", 
    name = marker_genes, 
    embedding = "UMAP",
    quantCut = c(0.01, 0.95),
    imputeWeights = getImputeWeights(proj)
)


plots <- lapply(p, function(x){
    x + guides(color = FALSE, fill = FALSE) + 
    theme_ArchR(baseSize = 6.5) +
    theme(plot.margin = unit(c(0, 0, 0, 0), "cm")) +
    theme(
        axis.text.x=element_blank(), 
        axis.ticks.x=element_blank(), 
        axis.text.y=element_blank(), 
        axis.ticks.y=element_blank()
    )
})
## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.

do.call(cowplot::plot_grid, c(list(ncol = 3),plots))

p <- plotEmbedding(
    ArchRProj = proj, 
    colorBy = "GeneExpressionMatrix", 
    name = "Gata6", 
    embedding = "UMAP",
    quantCut = c(0.01, 0.95)
)

p

Gene Expression Matrix UMAPs

rm(proj)
rm(new_scores)
## Warning in rm(new_scores): object 'new_scores' not found
gc(reset = TRUE)
##            used  (Mb) gc trigger  (Mb) max used  (Mb)
## Ncells  6654386 355.4   14774602 789.1  6654386 355.4
## Vcells 37397854 285.4  115950360 884.7 37397854 285.4
proj <- loadArchRProject("12_Ricards_peaks_ChromVar/")
scores <- getMatrixFromProject(proj, useMatrix = "GeneExpressionMatrix")
new_scores <- assays(scores)[[1]]
rownames(new_scores) <- rowData(scores)$name
rm(scores)
gc(reset = TRUE)
##             used   (Mb) gc trigger   (Mb)  max used   (Mb)
## Ncells   6984149  373.0   14774602  789.1   6984149  373.0
## Vcells 318918374 2433.2 1062600924 8107.1 318918374 2433.2

#marker_genes <- c("Hba-a2", "Hba-a1",  "Hbb-bh1")
p <- plotEmbedding(
    ArchRProj = proj, 
    colorBy = "GeneExpressionMatrix", 
    name = marker_genes, 
    embedding = "UMAP",
    quantCut = c(0.01, 0.95),
    imputeWeights = getImputeWeights(proj)
)


plots <- lapply(p, function(x){
    x + guides(color = FALSE, fill = FALSE) + 
    theme_ArchR(baseSize = 6.5) +
    theme(plot.margin = unit(c(0, 0, 0, 0), "cm")) +
    theme(
        axis.text.x=element_blank(), 
        axis.ticks.x=element_blank(), 
        axis.text.y=element_blank(), 
        axis.ticks.y=element_blank()
    )
})
## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.
## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.
## `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> = "none")` instead.

do.call(cowplot::plot_grid, c(list(ncol = 3),plots))

Correlation between gene activity score and gene expression:

Create cell aggregates

gene_expression <- getMatrixFromProject(proj, useMatrix = "GeneExpressionMatrix")
# extract gene expression matrix from gene expression object
gene_expr_mat <- assays(gene_expression)[[1]]
rownames(gene_expr_mat) <- rowData(gene_expression)$name
rm(gene_expression)
gc(reset = TRUE)
##             used   (Mb) gc trigger    (Mb)  max used   (Mb)
## Ncells   6987346  373.2   14774602   789.1   6987346  373.2
## Vcells 591591820 4513.5 1530321329 11675.5 591591820 4513.5

# randomly sample 5% of cells
sample_frac <- .05

# out of all cells in the gene expression matrix randomly sample 5% without replacement
sampled_cells <- sample(1:ncol(gene_expr_mat), sample_frac * ncol(gene_expr_mat), replace = FALSE)
print(paste0("Sampling ", length(sampled_cells), " cells."))
## [1] "Sampling 2299 cells."

Add PCA

library(Seurat)
#atac_seurat <- readRDS("atac_Seurat_object")
rna_seurat <- readRDS("Seurat_objects/rna_Seurat_object")

# get the cell embeddings
cell_embeddings <- Embeddings(object = rna_seurat, reduction = "pca")#[1:5, 1:5]


# get gene score matrix
#proj <- loadArchRProject("12_Ricards_peaks_ChromVar/")
#gene_scores <- getMatrixFromProject(proj, useMatrix = "GeneScoreMatrix")

First, we want to use only genes which we find in our newly computed gene activity scores from p2g links, the gene expression matrix and the imputed gene scores from scATAC.


#scores <- read.table("gene_scores_from_p2g/norm_simple_correction_no_scaling")
gc(reset = TRUE)
##              used    (Mb) gc trigger    (Mb)   max used    (Mb)
## Ncells    7656537   409.0   14774602   789.1    7656537   409.0
## Vcells 2648601521 20207.3 3800379050 28994.6 2648601521 20207.3
# create nearest neighbor graph
k <- 50

# use only cells we also have in gene expr/score matrix and the first 30 PCs
nn <- RANN::nn2(cell_embeddings[colnames(gene_expr_mat), 1:30], k = k)

# get genes of interest, the top 2000 most highly variable features
hvg <- rna_seurat@assays$originalexp@var.features

# we are left with 1772 genes
hvg <- hvg[hvg %in% rownames(new_scores)]
length(hvg)
## [1] 1772


# initialize matrix to store aggregate values
rna_score <- matrix(0, nrow = length(hvg), ncol = length(sampled_cells))
#atac_score <- matrix(0, nrow = length(hvg), ncol = length(sampled_cells))
p2g_score <- matrix(0, nrow = length(hvg), ncol = length(sampled_cells))
rownames(rna_score) <- hvg
#rownames(atac_score) <- hvg
rownames(p2g_score) <- hvg

# reorder columns to match gene expression matrix columns
#reord_gene_scores <- gene_score_mat[, colnames(gene_expr_mat)]
reord_p2g_scores <- new_scores[, colnames(gene_expr_mat)]
# 
# 
for(i in 1:length(sampled_cells)) {
# check progress
if (i %% 100 == 0) print(i)

# get current cell
cell <- sampled_cells[i]

# aggregate the rna score for one gene and all nearest neighbors of the first
# cell in our aggregate group of cells
rna_score[hvg, i] <- sparseMatrixStats::rowMeans2(
  gene_expr_mat[hvg, nn$nn.idx[cell, ]]
)

  # aggregate the atac score
p2g_score[hvg,i] <- sparseMatrixStats::rowMeans2(
  reord_p2g_scores[hvg, nn$nn.idx[cell,]]
  )
}
## [1] 100
## [1] 200
## [1] 300
## [1] 400
## [1] 500
## [1] 600
## [1] 700
## [1] 800
## [1] 900
## [1] 1000
## [1] 1100
## [1] 1200
## [1] 1300
## [1] 1400
## [1] 1500
## [1] 1600
## [1] 1700
## [1] 1800
## [1] 1900
## [1] 2000
## [1] 2100
## [1] 2200

# write.table(rna_score, "gene_scores_from_p2g/wrong_computation_good/rna_score")
# write.table(p2g_score, "gene_scores_from_p2g/wrong_computation_good/p2g_score")
rm(rna_seurat)
rm(gene_expr_mat)
gc(reset= TRUE)
##             used   (Mb) gc trigger    (Mb)  max used   (Mb)
## Ncells   7684871  410.5   14774602   789.1   7684871  410.5
## Vcells 606583433 4627.9 3040303240 23195.7 606583433 4627.9
hvg <- hvg[hvg %in%rownames(rna_score)]

comp_stats <- rbind(rna_score %>% set_rownames(paste0(hvg, "_#_rna_score")),
                    p2g_score %>% set_rownames(paste0(hvg, "_#_p2g_score"))) %>% 
  t() %>% 
  as.matrix() %>% 
  tibble::as_tibble() %>% 
  mutate(cell_name = sampled_cells) %>% 
  pivot_longer(cols = !cell_name, names_to = "name", values_to = "vals") %>%
  separate(name, sep = "_#_", into= c("gene", "statistic")) %>% 
  pivot_wider(names_from = statistic, values_from = vals)
marker_genes <- c("Lamb1",  "Sparc", "Elf5", "Ascl2", "Tfap2c", "Ttr",
                  "Apoa2", "Apoe", "Cystm1", "Emb", "Spink1",  "Krt19",
                  "Dkk1", "Grhl3", "Trp63", "Grhl2",  "Pax6", "Pax2",
                  "En1", "Foxd3", "Tfap2a", "Pax3", "Sox9",
                  "Six3", "Hesx1", "Irx3", "Sox2", "Hoxb9", "Cdx4", 
                  "Hes3", "Hba-a2", "Hba-a1",  "Hbb-bh1", "Gata1", "Cited4",
                   "Cdh5", "Pecam1", "Anxa5", "Etv2", "Igf2",
                  "Krt8", "Krt18", "Pmp22", "Ahnak", "Bmp4", "Tbx4", "Hoxa11", 
                  "Hoxa10", "Tnnt2", "Myl4",  "Myl7", "Acta2", 
                  "Smarcd3", "Tcf21", "Tbx6", "Dll1", "Aldh1a2", "Tcf15", 
                  "Meox1", "Tbx1", "Gbx2", "Cdx1", "Hoxb1", "Hes7", "Osr1", 
                  "Mesp2", "Lefty2", "Mesp1", "Cer1",  "Chrd", "T", 
                  "Foxa2", "Pax7", "Fgf8", "Lhx1", "Gsc", "Mixl1", "Otx2", "Hhex",
                   "Ifitm3", "Nkx1-2", "Eomes", "Nanog", "Utf1", 
                  "Epcam", "Pou5f1")


comp_stats %>%
  filter(gene %in% marker_genes) %>%
  ggplot(aes(x = rna_score, y = p2g_score)) +
  geom_point(size=0.4, alpha=0.4) +
  geom_smooth(formula = y ~ x, method = "lm", size=0.1) +
  #ggpubr::stat_cor(aes(label = ..r.label..), method = "pearson", r.digits=2) +
  facet_wrap(~ gene, scales="free", ncol = 5) +
  labs(x = "RNA Counts", y = "P2G Gene Activity Score")

LS0tCnRpdGxlOiAicDJnIHNjb3JlcyIKb3V0cHV0OiAKICBodG1sX2RvY3VtZW50OgogICAgdG9jOiB0cnVlCiAgICB0b2NfZGVwdGg6IDUKICAgIGNvZGVfZm9sZGluZzogaGlkZQogICAgdG9jX2Zsb2F0OiB0cnVlCiAgICBjb2RlX2Rvd25sb2FkOiB0cnVlCiAgICB0aGVtZTogY29zbW8KICAgIGhpZ2hsaWdodDogdGV4dG1hdGUKLS0tCgo8c3R5bGU+CmJvZHkgewp0ZXh0LWFsaWduOiBqdXN0aWZ5fQo8L3N0eWxlPgoKYGBge3Igc2V0dXAsIGluY2x1ZGU9RkFMU0V9CmtuaXRyOjpvcHRzX2NodW5rJHNldChjYWNoZSA9IEZBTFNFLCBhdXRvZGVwID0gVFJVRSwgCiAgICAgICAgICAgICAgICAgICAgICBjb2xsYXBzZSA9IFRSVUUsIG1lc3NhZ2UgPSBGQUxTRSkKa25pdHI6Om9wdHNfa25pdCRzZXQocm9vdC5kaXIgPSAiL29taWNzL2dyb3Vwcy9PRTA1MzMvaW50ZXJuYWwva2F0aGFyaW5hL3NjRG9SSS9nYXN0cnVsYXRpb25fZGF0YS8iKQpzZXR3ZCgiL29taWNzL2dyb3Vwcy9PRTA1MzMvaW50ZXJuYWwva2F0aGFyaW5hL3NjRG9SSS9nYXN0cnVsYXRpb25fZGF0YS8iKQoKc2V0LnNlZWQoMSkKYGBgCgogIApgYGB7cn0KIyMgTG9hZCBsaWJyYXJpZXMKCgpzdXBwcmVzc1BhY2thZ2VTdGFydHVwTWVzc2FnZXMoewpsaWJyYXJ5KEFyY2hSKQpsaWJyYXJ5KHJoZGY1KQpsaWJyYXJ5KHRpZHl2ZXJzZSkKbGlicmFyeShyZXRpY3VsYXRlKQpsaWJyYXJ5KHplbGxrb252ZXJ0ZXIpCmxpYnJhcnkoTWF0cml4KQpsaWJyYXJ5KGRpY2hyb21hdCkKI2xpYnJhcnkoY2FyZXQpCmg1ZGlzYWJsZUZpbGVMb2NraW5nKCl9KQpgYGAKCgoKYGBge3J9CnByb2ogPC0gbG9hZEFyY2hSUHJvamVjdCgiMTVfZ2VuZV9zY29yZXNfZnJvbV9wMmdfYXNfZ2VuZV9leHByX21hdHJpeC8iKQoKbWFya2VyX2dlbmVzIDwtIGMoIkxhbWIxIiwgICJTcGFyYyIsICJFbGY1IiwgIkFzY2wyIiwgIlRmYXAyYyIsICJUdHIiLAogICAgICAgICAgICAgICAgICAiQXBvYTIiLCAiQXBvZSIsICJDeXN0bTEiLCAiRW1iIiwgIlNwaW5rMSIsICAiS3J0MTkiLAogICAgICAgICAgICAgICAgICAiRGtrMSIsICJHcmhsMyIsICJUcnA2MyIsICJHcmhsMiIsICAiUGF4NiIsICJQYXgyIiwKICAgICAgICAgICAgICAgICAgIkVuMSIsICJGb3hkMyIsICJUZmFwMmEiLCAiUGF4MyIsICJTb3g5IiwKICAgICAgICAgICAgICAgICAgIlNpeDMiLCAiSGVzeDEiLCAiSXJ4MyIsICJTb3gyIiwgIkhveGI5IiwgIkNkeDQiLCAKICAgICAgICAgICAgICAgICAgIkhlczMiLCAiSGJhLWEyIiwgIkhiYS1hMSIsICAiSGJiLWJoMSIsICJHYXRhMSIsICJDaXRlZDQiLAogICAgICAgICAgICAgICAgICAgIkNkaDUiLCAiUGVjYW0xIiwgIkFueGE1IiwgIkV0djIiLCAiSWdmMiIsCiAgICAgICAgICAgICAgICAgICJLcnQ4IiwgIktydDE4IiwgIlBtcDIyIiwgIkFobmFrIiwgIkJtcDQiLCAiVGJ4NCIsICJIb3hhMTEiLCAKICAgICAgICAgICAgICAgICAgIkhveGExMCIsICJUbm50MiIsICJNeWw0IiwgICJNeWw3IiwgIkFjdGEyIiwgCiAgICAgICAgICAgICAgICAgICJTbWFyY2QzIiwgIlRjZjIxIiwgIlRieDYiLCAiRGxsMSIsICJBbGRoMWEyIiwgIlRjZjE1IiwgCiAgICAgICAgICAgICAgICAgICJNZW94MSIsICJUYngxIiwgIkdieDIiLCAiQ2R4MSIsICJIb3hiMSIsICJIZXM3IiwgIk9zcjEiLCAKICAgICAgICAgICAgICAgICAgIk1lc3AyIiwgIkxlZnR5MiIsICJNZXNwMSIsICJDZXIxIiwgICJDaHJkIiwgIlQiLCAKICAgICAgICAgICAgICAgICAgIkZveGEyIiwgIlBheDciLCAiRmdmOCIsICJMaHgxIiwgIkdzYyIsICJNaXhsMSIsICJPdHgyIiwgIkhoZXgiLAogICAgICAgICAgICAgICAgICAgIklmaXRtMyIsICJOa3gxLTIiLCAiRW9tZXMiLCAiTmFub2ciLCAiVXRmMSIsIAogICAgICAgICAgICAgICAgICAiRXBjYW0iLCAiUG91NWYxIikKYGBgCgojIyMgR2VuZSBTY29yZXMgVU1BUHMKCmBgYHtyLCBmaWcud2lkdGg9MTAsIGZpZy5oZWlnaHQ9NTB9CgojbWFya2VyX2dlbmVzIDwtIGMoIkhiYS1hMiIsICJIYmEtYTEiLCAgIkhiYi1iaDEiKQpwIDwtIHBsb3RFbWJlZGRpbmcoCiAgICBBcmNoUlByb2ogPSBwcm9qLCAKICAgIGNvbG9yQnkgPSAiR2VuZUV4cHJlc3Npb25NYXRyaXgiLCAKICAgIG5hbWUgPSBtYXJrZXJfZ2VuZXMsIAogICAgZW1iZWRkaW5nID0gIlVNQVAiLAogICAgcXVhbnRDdXQgPSBjKDAuMDEsIDAuOTUpLAogICAgaW1wdXRlV2VpZ2h0cyA9IGdldEltcHV0ZVdlaWdodHMocHJvaikKKQoKCnBsb3RzIDwtIGxhcHBseShwLCBmdW5jdGlvbih4KXsKICAgIHggKyBndWlkZXMoY29sb3IgPSBGQUxTRSwgZmlsbCA9IEZBTFNFKSArIAogICAgdGhlbWVfQXJjaFIoYmFzZVNpemUgPSA2LjUpICsKICAgIHRoZW1lKHBsb3QubWFyZ2luID0gdW5pdChjKDAsIDAsIDAsIDApLCAiY20iKSkgKwogICAgdGhlbWUoCiAgICAgICAgYXhpcy50ZXh0Lng9ZWxlbWVudF9ibGFuaygpLCAKICAgICAgICBheGlzLnRpY2tzLng9ZWxlbWVudF9ibGFuaygpLCAKICAgICAgICBheGlzLnRleHQueT1lbGVtZW50X2JsYW5rKCksIAogICAgICAgIGF4aXMudGlja3MueT1lbGVtZW50X2JsYW5rKCkKICAgICkKfSkKCmRvLmNhbGwoY293cGxvdDo6cGxvdF9ncmlkLCBjKGxpc3QobmNvbCA9IDMpLHBsb3RzKSkKCmBgYAoKYGBgI3tyfQpwIDwtIHBsb3RFbWJlZGRpbmcoCiAgICBBcmNoUlByb2ogPSBwcm9qLCAKICAgIGNvbG9yQnkgPSAiR2VuZUV4cHJlc3Npb25NYXRyaXgiLCAKICAgIG5hbWUgPSAiR2F0YTYiLCAKICAgIGVtYmVkZGluZyA9ICJVTUFQIiwKICAgIHF1YW50Q3V0ID0gYygwLjAxLCAwLjk1KQopCgpwCmBgYAoKIyMjIEdlbmUgRXhwcmVzc2lvbiBNYXRyaXggVU1BUHMKYGBge3J9CnJtKHByb2opCnJtKG5ld19zY29yZXMpCmdjKHJlc2V0ID0gVFJVRSkKcHJvaiA8LSBsb2FkQXJjaFJQcm9qZWN0KCIxMl9SaWNhcmRzX3BlYWtzX0Nocm9tVmFyLyIpCnNjb3JlcyA8LSBnZXRNYXRyaXhGcm9tUHJvamVjdChwcm9qLCB1c2VNYXRyaXggPSAiR2VuZUV4cHJlc3Npb25NYXRyaXgiKQpuZXdfc2NvcmVzIDwtIGFzc2F5cyhzY29yZXMpW1sxXV0Kcm93bmFtZXMobmV3X3Njb3JlcykgPC0gcm93RGF0YShzY29yZXMpJG5hbWUKcm0oc2NvcmVzKQpnYyhyZXNldCA9IFRSVUUpCgojbWFya2VyX2dlbmVzIDwtIGMoIkhiYS1hMiIsICJIYmEtYTEiLCAgIkhiYi1iaDEiKQpwIDwtIHBsb3RFbWJlZGRpbmcoCiAgICBBcmNoUlByb2ogPSBwcm9qLCAKICAgIGNvbG9yQnkgPSAiR2VuZUV4cHJlc3Npb25NYXRyaXgiLCAKICAgIG5hbWUgPSBtYXJrZXJfZ2VuZXMsIAogICAgZW1iZWRkaW5nID0gIlVNQVAiLAogICAgcXVhbnRDdXQgPSBjKDAuMDEsIDAuOTUpLAogICAgaW1wdXRlV2VpZ2h0cyA9IGdldEltcHV0ZVdlaWdodHMocHJvaikKKQoKCnBsb3RzIDwtIGxhcHBseShwLCBmdW5jdGlvbih4KXsKICAgIHggKyBndWlkZXMoY29sb3IgPSBGQUxTRSwgZmlsbCA9IEZBTFNFKSArIAogICAgdGhlbWVfQXJjaFIoYmFzZVNpemUgPSA2LjUpICsKICAgIHRoZW1lKHBsb3QubWFyZ2luID0gdW5pdChjKDAsIDAsIDAsIDApLCAiY20iKSkgKwogICAgdGhlbWUoCiAgICAgICAgYXhpcy50ZXh0Lng9ZWxlbWVudF9ibGFuaygpLCAKICAgICAgICBheGlzLnRpY2tzLng9ZWxlbWVudF9ibGFuaygpLCAKICAgICAgICBheGlzLnRleHQueT1lbGVtZW50X2JsYW5rKCksIAogICAgICAgIGF4aXMudGlja3MueT1lbGVtZW50X2JsYW5rKCkKICAgICkKfSkKCmRvLmNhbGwoY293cGxvdDo6cGxvdF9ncmlkLCBjKGxpc3QobmNvbCA9IDMpLHBsb3RzKSkKCgpgYGAKCgoKIyBDb3JyZWxhdGlvbiBiZXR3ZWVuIGdlbmUgYWN0aXZpdHkgc2NvcmUgYW5kIGdlbmUgZXhwcmVzc2lvbjoKCiMjIENyZWF0ZSBjZWxsIGFnZ3JlZ2F0ZXMKCmBgYHtyfQpnZW5lX2V4cHJlc3Npb24gPC0gZ2V0TWF0cml4RnJvbVByb2plY3QocHJvaiwgdXNlTWF0cml4ID0gIkdlbmVFeHByZXNzaW9uTWF0cml4IikKIyBleHRyYWN0IGdlbmUgZXhwcmVzc2lvbiBtYXRyaXggZnJvbSBnZW5lIGV4cHJlc3Npb24gb2JqZWN0CmdlbmVfZXhwcl9tYXQgPC0gYXNzYXlzKGdlbmVfZXhwcmVzc2lvbilbWzFdXQpyb3duYW1lcyhnZW5lX2V4cHJfbWF0KSA8LSByb3dEYXRhKGdlbmVfZXhwcmVzc2lvbikkbmFtZQpybShnZW5lX2V4cHJlc3Npb24pCmdjKHJlc2V0ID0gVFJVRSkKCiMgcmFuZG9tbHkgc2FtcGxlIDUlIG9mIGNlbGxzCnNhbXBsZV9mcmFjIDwtIC4wNQoKIyBvdXQgb2YgYWxsIGNlbGxzIGluIHRoZSBnZW5lIGV4cHJlc3Npb24gbWF0cml4IHJhbmRvbWx5IHNhbXBsZSA1JSB3aXRob3V0IHJlcGxhY2VtZW50CnNhbXBsZWRfY2VsbHMgPC0gc2FtcGxlKDE6bmNvbChnZW5lX2V4cHJfbWF0KSwgc2FtcGxlX2ZyYWMgKiBuY29sKGdlbmVfZXhwcl9tYXQpLCByZXBsYWNlID0gRkFMU0UpCnByaW50KHBhc3RlMCgiU2FtcGxpbmcgIiwgbGVuZ3RoKHNhbXBsZWRfY2VsbHMpLCAiIGNlbGxzLiIpKQpgYGAKCgoKIyMgQWRkIFBDQQoKCgoKYGBge3J9CmxpYnJhcnkoU2V1cmF0KQojYXRhY19zZXVyYXQgPC0gcmVhZFJEUygiYXRhY19TZXVyYXRfb2JqZWN0IikKcm5hX3NldXJhdCA8LSByZWFkUkRTKCJTZXVyYXRfb2JqZWN0cy9ybmFfU2V1cmF0X29iamVjdCIpCgojIGdldCB0aGUgY2VsbCBlbWJlZGRpbmdzCmNlbGxfZW1iZWRkaW5ncyA8LSBFbWJlZGRpbmdzKG9iamVjdCA9IHJuYV9zZXVyYXQsIHJlZHVjdGlvbiA9ICJwY2EiKSNbMTo1LCAxOjVdCgoKIyBnZXQgZ2VuZSBzY29yZSBtYXRyaXgKI3Byb2ogPC0gbG9hZEFyY2hSUHJvamVjdCgiMTJfUmljYXJkc19wZWFrc19DaHJvbVZhci8iKQojZ2VuZV9zY29yZXMgPC0gZ2V0TWF0cml4RnJvbVByb2plY3QocHJvaiwgdXNlTWF0cml4ID0gIkdlbmVTY29yZU1hdHJpeCIpCmBgYAoKCkZpcnN0LCB3ZSB3YW50IHRvIHVzZSBvbmx5IGdlbmVzIHdoaWNoIHdlIGZpbmQgaW4gb3VyIG5ld2x5IGNvbXB1dGVkIGdlbmUKYWN0aXZpdHkgc2NvcmVzIGZyb20gcDJnIGxpbmtzLCB0aGUgZ2VuZSBleHByZXNzaW9uIG1hdHJpeCBhbmQgdGhlIGltcHV0ZWQgZ2VuZQpzY29yZXMgZnJvbSBzY0FUQUMuCgpgYGB7cn0KCiNzY29yZXMgPC0gcmVhZC50YWJsZSgiZ2VuZV9zY29yZXNfZnJvbV9wMmcvbm9ybV9zaW1wbGVfY29ycmVjdGlvbl9ub19zY2FsaW5nIikKYGBgCgoKYGBge3J9CmdjKHJlc2V0ID0gVFJVRSkKIyBjcmVhdGUgbmVhcmVzdCBuZWlnaGJvciBncmFwaAprIDwtIDUwCgojIHVzZSBvbmx5IGNlbGxzIHdlIGFsc28gaGF2ZSBpbiBnZW5lIGV4cHIvc2NvcmUgbWF0cml4IGFuZCB0aGUgZmlyc3QgMzAgUENzCm5uIDwtIFJBTk46Om5uMihjZWxsX2VtYmVkZGluZ3NbY29sbmFtZXMoZ2VuZV9leHByX21hdCksIDE6MzBdLCBrID0gaykKCiMgZ2V0IGdlbmVzIG9mIGludGVyZXN0LCB0aGUgdG9wIDIwMDAgbW9zdCBoaWdobHkgdmFyaWFibGUgZmVhdHVyZXMKaHZnIDwtIHJuYV9zZXVyYXRAYXNzYXlzJG9yaWdpbmFsZXhwQHZhci5mZWF0dXJlcwoKIyB3ZSBhcmUgbGVmdCB3aXRoIDE3NzIgZ2VuZXMKaHZnIDwtIGh2Z1todmcgJWluJSByb3duYW1lcyhuZXdfc2NvcmVzKV0KbGVuZ3RoKGh2ZykKCgojIGluaXRpYWxpemUgbWF0cml4IHRvIHN0b3JlIGFnZ3JlZ2F0ZSB2YWx1ZXMKcm5hX3Njb3JlIDwtIG1hdHJpeCgwLCBucm93ID0gbGVuZ3RoKGh2ZyksIG5jb2wgPSBsZW5ndGgoc2FtcGxlZF9jZWxscykpCiNhdGFjX3Njb3JlIDwtIG1hdHJpeCgwLCBucm93ID0gbGVuZ3RoKGh2ZyksIG5jb2wgPSBsZW5ndGgoc2FtcGxlZF9jZWxscykpCnAyZ19zY29yZSA8LSBtYXRyaXgoMCwgbnJvdyA9IGxlbmd0aChodmcpLCBuY29sID0gbGVuZ3RoKHNhbXBsZWRfY2VsbHMpKQpyb3duYW1lcyhybmFfc2NvcmUpIDwtIGh2Zwojcm93bmFtZXMoYXRhY19zY29yZSkgPC0gaHZnCnJvd25hbWVzKHAyZ19zY29yZSkgPC0gaHZnCgojIHJlb3JkZXIgY29sdW1ucyB0byBtYXRjaCBnZW5lIGV4cHJlc3Npb24gbWF0cml4IGNvbHVtbnMKI3Jlb3JkX2dlbmVfc2NvcmVzIDwtIGdlbmVfc2NvcmVfbWF0WywgY29sbmFtZXMoZ2VuZV9leHByX21hdCldCnJlb3JkX3AyZ19zY29yZXMgPC0gbmV3X3Njb3Jlc1ssIGNvbG5hbWVzKGdlbmVfZXhwcl9tYXQpXQojIAojIApmb3IoaSBpbiAxOmxlbmd0aChzYW1wbGVkX2NlbGxzKSkgewojIGNoZWNrIHByb2dyZXNzCmlmIChpICUlIDEwMCA9PSAwKSBwcmludChpKQoKIyBnZXQgY3VycmVudCBjZWxsCmNlbGwgPC0gc2FtcGxlZF9jZWxsc1tpXQoKIyBhZ2dyZWdhdGUgdGhlIHJuYSBzY29yZSBmb3Igb25lIGdlbmUgYW5kIGFsbCBuZWFyZXN0IG5laWdoYm9ycyBvZiB0aGUgZmlyc3QKIyBjZWxsIGluIG91ciBhZ2dyZWdhdGUgZ3JvdXAgb2YgY2VsbHMKcm5hX3Njb3JlW2h2ZywgaV0gPC0gc3BhcnNlTWF0cml4U3RhdHM6OnJvd01lYW5zMigKICBnZW5lX2V4cHJfbWF0W2h2Zywgbm4kbm4uaWR4W2NlbGwsIF1dCikKCiAgIyBhZ2dyZWdhdGUgdGhlIGF0YWMgc2NvcmUKcDJnX3Njb3JlW2h2ZyxpXSA8LSBzcGFyc2VNYXRyaXhTdGF0czo6cm93TWVhbnMyKAogIHJlb3JkX3AyZ19zY29yZXNbaHZnLCBubiRubi5pZHhbY2VsbCxdXQogICkKfQoKIyB3cml0ZS50YWJsZShybmFfc2NvcmUsICJnZW5lX3Njb3Jlc19mcm9tX3AyZy93cm9uZ19jb21wdXRhdGlvbl9nb29kL3JuYV9zY29yZSIpCiMgd3JpdGUudGFibGUocDJnX3Njb3JlLCAiZ2VuZV9zY29yZXNfZnJvbV9wMmcvd3JvbmdfY29tcHV0YXRpb25fZ29vZC9wMmdfc2NvcmUiKQoKYGBgCgoKYGBge3J9CnJtKHJuYV9zZXVyYXQpCnJtKGdlbmVfZXhwcl9tYXQpCmdjKHJlc2V0PSBUUlVFKQpgYGAKCgpgYGB7cn0KaHZnIDwtIGh2Z1todmcgJWluJXJvd25hbWVzKHJuYV9zY29yZSldCgpjb21wX3N0YXRzIDwtIHJiaW5kKHJuYV9zY29yZSAlPiUgc2V0X3Jvd25hbWVzKHBhc3RlMChodmcsICJfI19ybmFfc2NvcmUiKSksCiAgICAgICAgICAgICAgICAgICAgcDJnX3Njb3JlICU+JSBzZXRfcm93bmFtZXMocGFzdGUwKGh2ZywgIl8jX3AyZ19zY29yZSIpKSkgJT4lIAogIHQoKSAlPiUgCiAgYXMubWF0cml4KCkgJT4lIAogIHRpYmJsZTo6YXNfdGliYmxlKCkgJT4lIAogIG11dGF0ZShjZWxsX25hbWUgPSBzYW1wbGVkX2NlbGxzKSAlPiUgCiAgcGl2b3RfbG9uZ2VyKGNvbHMgPSAhY2VsbF9uYW1lLCBuYW1lc190byA9ICJuYW1lIiwgdmFsdWVzX3RvID0gInZhbHMiKSAlPiUKICBzZXBhcmF0ZShuYW1lLCBzZXAgPSAiXyNfIiwgaW50bz0gYygiZ2VuZSIsICJzdGF0aXN0aWMiKSkgJT4lIAogIHBpdm90X3dpZGVyKG5hbWVzX2Zyb20gPSBzdGF0aXN0aWMsIHZhbHVlc19mcm9tID0gdmFscykKYGBgCgoKYGBge3IsIGZpZy53aWRodCA9IDUsIGZpZy5oZWlnaHQgPSAyMH0KbWFya2VyX2dlbmVzIDwtIGMoIkxhbWIxIiwgICJTcGFyYyIsICJFbGY1IiwgIkFzY2wyIiwgIlRmYXAyYyIsICJUdHIiLAogICAgICAgICAgICAgICAgICAiQXBvYTIiLCAiQXBvZSIsICJDeXN0bTEiLCAiRW1iIiwgIlNwaW5rMSIsICAiS3J0MTkiLAogICAgICAgICAgICAgICAgICAiRGtrMSIsICJHcmhsMyIsICJUcnA2MyIsICJHcmhsMiIsICAiUGF4NiIsICJQYXgyIiwKICAgICAgICAgICAgICAgICAgIkVuMSIsICJGb3hkMyIsICJUZmFwMmEiLCAiUGF4MyIsICJTb3g5IiwKICAgICAgICAgICAgICAgICAgIlNpeDMiLCAiSGVzeDEiLCAiSXJ4MyIsICJTb3gyIiwgIkhveGI5IiwgIkNkeDQiLCAKICAgICAgICAgICAgICAgICAgIkhlczMiLCAiSGJhLWEyIiwgIkhiYS1hMSIsICAiSGJiLWJoMSIsICJHYXRhMSIsICJDaXRlZDQiLAogICAgICAgICAgICAgICAgICAgIkNkaDUiLCAiUGVjYW0xIiwgIkFueGE1IiwgIkV0djIiLCAiSWdmMiIsCiAgICAgICAgICAgICAgICAgICJLcnQ4IiwgIktydDE4IiwgIlBtcDIyIiwgIkFobmFrIiwgIkJtcDQiLCAiVGJ4NCIsICJIb3hhMTEiLCAKICAgICAgICAgICAgICAgICAgIkhveGExMCIsICJUbm50MiIsICJNeWw0IiwgICJNeWw3IiwgIkFjdGEyIiwgCiAgICAgICAgICAgICAgICAgICJTbWFyY2QzIiwgIlRjZjIxIiwgIlRieDYiLCAiRGxsMSIsICJBbGRoMWEyIiwgIlRjZjE1IiwgCiAgICAgICAgICAgICAgICAgICJNZW94MSIsICJUYngxIiwgIkdieDIiLCAiQ2R4MSIsICJIb3hiMSIsICJIZXM3IiwgIk9zcjEiLCAKICAgICAgICAgICAgICAgICAgIk1lc3AyIiwgIkxlZnR5MiIsICJNZXNwMSIsICJDZXIxIiwgICJDaHJkIiwgIlQiLCAKICAgICAgICAgICAgICAgICAgIkZveGEyIiwgIlBheDciLCAiRmdmOCIsICJMaHgxIiwgIkdzYyIsICJNaXhsMSIsICJPdHgyIiwgIkhoZXgiLAogICAgICAgICAgICAgICAgICAgIklmaXRtMyIsICJOa3gxLTIiLCAiRW9tZXMiLCAiTmFub2ciLCAiVXRmMSIsIAogICAgICAgICAgICAgICAgICAiRXBjYW0iLCAiUG91NWYxIikKCgpjb21wX3N0YXRzICU+JQogIGZpbHRlcihnZW5lICVpbiUgbWFya2VyX2dlbmVzKSAlPiUKICBnZ3Bsb3QoYWVzKHggPSBybmFfc2NvcmUsIHkgPSBwMmdfc2NvcmUpKSArCiAgZ2VvbV9wb2ludChzaXplPTAuNCwgYWxwaGE9MC40KSArCiAgZ2VvbV9zbW9vdGgoZm9ybXVsYSA9IHkgfiB4LCBtZXRob2QgPSAibG0iLCBzaXplPTAuMSkgKwogICNnZ3B1YnI6OnN0YXRfY29yKGFlcyhsYWJlbCA9IC4uci5sYWJlbC4uKSwgbWV0aG9kID0gInBlYXJzb24iLCByLmRpZ2l0cz0yKSArCiAgZmFjZXRfd3JhcCh+IGdlbmUsIHNjYWxlcz0iZnJlZSIsIG5jb2wgPSA1KSArCiAgbGFicyh4ID0gIlJOQSBDb3VudHMiLCB5ID0gIlAyRyBHZW5lIEFjdGl2aXR5IFNjb3JlIikKYGBgCg==